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ABSTRACT 

We present theoretical calculations for the differential distribution of stellar orbital 
eccentricity in a galaxy halo, assuming that the stars constitute a spherical, collisionless 
system in dynamical equilibrium with a dark matter halo. In order to define the 
eccentricity e of a halo star for given energy E and angular momentum L, we adopt 
two types of gravitational potential, such as an isochrone potential and a Navarro- 
Frenk- White potential, that could form two ends covering in-between any realistic 
potential of dark matter halo. Based on a distribution function of the form f(E, L) 
that allows constant anisotropy in velocity dispersions characterized by a parameter 
j3, we find that the eccentricity distribution is a monotonically increasing function of e 
for the case of highly radially anisotropic velocity dispersions (f3 > 0.6), while showing 
a hump-like shape for the cases from radial through tangential velocity anisotropy 
(j3 < 0.6). We also find that when the velocity anisotropy agrees with that observed 
for the Milky Way halo stars (j3 ~ 0.5 — 0.7), a nearly linear eccentricity distribution 
of N(e) oc e results at e < 0.7, largely independent of the potential adopted. Our 
theoretical eccentricity distribution would be a vital tool of examining how far out 
in the halo the dynamical equilibrium has been achieved, through comparison with 
kinematics of halo stars sampled at greater distances. Given that large surveys of the 
SEGUE and Gaia projects would be in progress, we discuss how our results would 
serve as a new guide in exploring the formation and evolution of the Milky Way halo. 

Key words: Galaxy: halo - Galaxy: kinematics and dynamics - Galaxy: formation 
- Galaxy: evolution - stellar dynamics - methods: analytical. 



1 INTRODUCTION 

Studies of large-scale structures in the universe and fluc- 
tuations in the cosmic microwave background strongly fa- 
vor a A-cold dark matter (ACDM) cosmology (e.g., Cole et 
al. 2005; Dunkley et al. 2009). The formation of structures 
in this cosmology is a process of hierarchical clustering, in 
the sense that numerous CDM lumps cluster gravitationally 
and merge together to form larger structures (White & Rees 
1978; Blumenthal et al. 1984). 

Dark halos of galaxy systems are similarly formed via 
clustering of subhalos as a result of CDM agglomerations 
that reach the maximum expansion then turn around to col- 
lapse in the background expanding medium, but a detailed 
process leading to the halo formation from primordial den- 
sity fluctuations is highly nonlinear and is not as simple as 
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the formation of larger structures in the universe (e.g., for 
review see Ostriker 1993 and Bertschinger 1998). 

High-resolution ACDM simulations for the halo forma- 
tion generically show that mergers and collisions of subhalos 
induce the overall collapse and virialize the inner region of 
host halo, while surviving subhalos orbit as separate enti- 
ties within the inner virialized region of halo (e.g., Moore et 
al. 1999; Ghigna et al. 2000; Helmi, White & Springel 2003; 
Valluri et al. 2007). A majority of stars formed through this 
build-up of halo are expected to have also experienced the re- 
distribution of energy and momentum that drives the phase 
mixing or violent relaxation towards the dynamical equilib- 
rium (Lynden-Bell 1967). This leads to an idea that a stellar 
halo, which can be regarded as a collisionless system, holds 
the dynamical information just after the last violent relax- 
ation in forming the halo. 

We then take an approach to find out the relics of the 
formation of the Milky Way halo from the kinematics of halo 
stars. Among many of their kinematic properties available at 
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present and in the near future, the differential distribution 
of stellar orbital eccentricity N(e) seems to be of special 
importance. The orbital eccentricity of a star is a quasi- 
adiabatic invariant (Eggen, Lynden-Bell & Sandage 1962; 
Lynden-Bell 1963) and is unaffected by the small and slow 
variation of the gravitational potential that might have oc- 
curred after the major formation of halo stars. It is therefore 
most likely that the shape of N(e) has been conserved un- 
til present. With this consideration, comparing the observed 
shape of iV(e) for halo stars with the theoretical one for the 
halo in dynamical equilibrium, we could explore how far out 
in the halo the dynamical equilibrium was achieved. Conse- 
quently, N(e) serves as a new test of halo formation scenario 
in a ACDM cosmology. 

As a useful way to derive N(e) theoretically, we con- 
sider the orbit of halo stars in assumed gravitational poten- 
tials of the halo. In section 2, we present our formulation 
to calculate N(e) under some plausible assumptions for the 
halo, and apply it to two extreme gravitational potentials of 
academic interest. The results for realistic cases are shown 
for the isochrone potential and for the Navarro-Frenk- White 
(NFW) potential in section 3. We summarize the results and 
discuss the prospects of investigating the formation and evo- 
lution of the Milky Way halo in section 4. 



2 FORMULATION 

We assume that the halo stars constitute a spherical, colli- 
sionless system in dynamical equilibrium with a dark halo. 
Since the dark matter is known to dominate the total mass 
of the galaxy system, the motion of halo stars is governed 
by the gravitational potential of dark halo. 

2.1 Stellar orbital eccentricity in a model halo 

When a spherical halo potential V(r) is given with respect to 
the galaxy center, the energy E and the angular momentum 
L of a star at the position r with the velocity v are written 
respectively as 

E=^v 2 + V{r), and L = \L\ = \r X v\ , (1) 

where r = \r\. The orbital eccentricity of a star is practically 
defined as 



< "r ^*pe 



(2) 



where r apo and r por i are the apo- and peri-centric distances, 
respectively, and are given by two real solutions (r apo > 
»~pcri) of the following equation: 

L 2 



E = V{r) + 



2r 2 



V eg (L;r). 



(3) 



It is evident from equations © and ((3} that a pair of (E, L) 
has a one-to-one correspondence to (E, e), but there is a re- 
gion of (E, L) in which two real solutions are not allowed and 
thus, except for the case of circular orbits, the eccentricity 
cannot be defined. Since such unbound orbits do not form a 
steady population of stellar halo, we neglect them and ex- 
clusively consider stars with bound orbits. Constraints on E 
and L that allow bound orbits are presented in Appendix 
A. 



2.2 Differential distribution of stellar orbital 
eccentricity 

Let f(r, v) be the distribution function of halo stars, then 
the number of halo stars in a phase space volume d 3 rd 3 v 
centered at (r,v) is given by f(r 1 v)d 3 rd 3 v. According to 
the strong Jeans theorem, the distribution function should 
be expressed in terms of isolating integrals only (Lynden- 
Bell 1960, 1962). For a spherical system that is invariant 
under rotation, it takes a form of either f(E) or f(E,L), 
depending on whether the stellar velocity dispersion is 
isotropic or anisotropic, respectively. 

The velocity dispersion observed for halo stars is radi- 
ally anisotropic (e.g., Yoshii & Saio 1979; Gilmore, Wyse, 
& Kuijken 1989). Furthermore, recent observations for halo 
stars within the distance of 10 kpc away from us show that 
the shape of velocity ellipsoid is constant and its principal 
axes are well aligned with the spherical coordinates (Carollo 
et al. 2007; Bond et al. 2009). If we extrapolate this fact to 
a whole system, one simple form of the distribution function 
is 



f(E,L) 



g{E)L- 211 , if (E,L) is 'allowed' 
0, otherwise, 



(4) 



where g(E) is a function of E (e.g., Binney & Tremine 2008). 
Here, j3 is a constant value of velocity anisotropy parameter 
defined as 



^1-i 



(5) 



where oy is the radial velocity dispersion and a t is the tan- 
gential velocity dispersion projected onto the spherical 6 — cj> 
surface. Although j3 is about 0.5 — 0.7 observationally (e.g., 
Bond et al. 2009; Smith et al. 2009; Carollo et al. 2010), we 
will use it as a constant parameter below. 

By changing the variables and integrating over the 
spherical coordinates, the number of stars in d 3 rd 3 v reduces 
to 



N{E,L)dEdL 2 = 4-k 2 g(E)L- 2l3 TrdEdL 2 , 
with the radial period of stellar orbit given by 

dr 



T r (E,L) 



dr 



,-, sj2[E-V ctt {L-r)\ 



(6) 



(7) 



Since L 2 is a function of E and e, we here introduce the 
^-dependent differential eccentricity distribution as 



np(E,e) =4nL-^ ■ T, 



de 



(8) 



We then express the differential eccentricity distribution as 



N (e) 



g(E)n (E,e)dE. 



(9) 



It is apparent from this equation that Np(e) is a weighted 
sum of np(E, e) with a weight function of g(E). Thus, once 
the gravitational potential V(r) and the velocity anisotropy 
parameter /3 are specified, we can formally obtain np(E,e), 
and also Np(e) after integrating np{E,e) over E with its 
appropriate weight. 
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2.3 Extreme cases of mass distribution 

In this subsection, mostly for pedagogical purpose, we con- 
sider two extreme cases of mass distribution such as the 
point mass at the center and the homogeneous distribution 
in the truncated sphere. These cases allow analytic expres- 
sion of n/3 (E,e), and because it is separable in E and e, 
Np(e) can also be obtained except for its normalization. 
Therefore, these cases are helpful to understand the results 
for any more realistic cases. 



2.3.1 Central point mass 

The gravitational potential arising from the central point 
mass is Keplerian and is given by 



V(r) = 



GM 



(10) 



where M is the total mass of dark halo and G is the gravita- 
tional constant. For bound orbits with E < 0, there are two 
real and positive solutions for equation ©, or equivalently, 



(-2E)r 2 - 2GMr + L 2 = 0. 



(11) 



The orbital eccentricity is expressed in terms of (E, L) as 



(-2E)L 2 



(GM) 2 ' 
and the other relevant quantities are neatly expressed as 



(12) 



T r = 



2-kGM , , 2 {GM) 2 M 



(-2E)i 



and L = 



-2E 



1-e 



(13) 



Substitution of these quantities in equations ((8| and ((9} 
gives the independent differential eccentricity distribution 



np(E,e) = l&-K A {GMy-' p {-2E) 



\p-i. 



(1-e*)"' 

and the differential eccentricity distribution 



(14) 



Np(e) 



16tt 3 (GM) 3 " 



/ 



g(E)(-2Ef-2dE 



(1 



,2)0 * 

(15) 



Since Np(e) oc np(E,e), we normalize Np(e) such that 
f Np(e)de = 1, and write 



normalized Np(e) = 2(1 - f}) 



(1-e 2 ) 



2\S ■ 



08^1). 



(16) 



The results of Np (e) for several values of /3 are shown on the 
left panel of Figure [T] For the case of j3 = (isotropic veloc- 
ity dispersion), Np(e) is exactly proportional to e (Binney & 
Tremaine 2008) and we call it the linear eccentricity distri- 
bution. For < /3 < 1 (radially anisotropic velocity disper- 
sion) , Np (e) is a rapidly increasing function of e with a peak 
always at e = 1. On the other hand, for f3 < (tangentially 
anisotropic velocity dispersion), Np(e) shows a hump-like 
e-distribution around a single peak at e = (1 — 2f3)~ 1 ' 2 . 

We should notice that the linear trend of Np oc e pre- 
vails in a range of < e < 0.3 regardless of /3, while the be- 
havior of Np (e) is very sensitive to j3 in a range of 0.6 < e < 1 
and the difference there clearly shows up. 



2.3.2 Truncated homogeneous sphere 

A homogeneous density distribution within truncated sphere 
is expressed as 



p(r) = 



3 A/ 
A'Kr'r 



0, 



if r < rt 

otherwise, 



(17) 



where M is the total mass of dark halo and r t is the trun- 
cation radius. The gravitational potential arising from this 
density distribution is given by 



V(r) = 



3GM , GM 
2n ~ l ~ 2n 



(*)' 



if r < r t 
otherwise. 



(18) 



We consider only stars with E < Et = —GM/rt, which 
guarantees the stars to be confined inside the truncated ra- 
dius rt. Thus, bound orbits within the truncated sphere are 
allowed if -E m in < E < E t where we note B m i n = (3/2) Et. 
In this limited range of E, there are two real and positive 
solutions for equation ©, or equivalently, 



GMrJ^A -2(E-E min )r 2 f^\ + L 2 = 0, 

if and only if 
0< D< 1, 
where 

GML 2 



rf(E-E min ) 2 ' 
The orbital eccentricity is expressed in terms of D as 



i-Vd 
T+7d' 



(19) 

(20) 
(21) 

(22) 



and the other relevant quantities are expressed in terms of 

(E, L) as 



T r 



GM' 



and L 



7t^ E - E " 



„2\ 2 



l + e s 



(23) 



Consequently, we obtain 

„3 \ 1-/9 



np(E,e) = 32tt' 



\GM) 



(E - E n 



-2/3 e(l 



„2n1-2/3 



(24) 



and 



Np(e) = 32tt : 



\gm) 



g(E)(E-E mln 

. e (l 



.3-3/8 



dE 



„2>l-2/3 



(1 + e 2 ) 3 " 2 ' 3 



(25) 



As in the point mass model, since Np{e) oc np(E,e), we 
normalize Np(e) such that J Np(e)de — 1, and write 

normalized Np{e) = 8(1 - /})-£ -^- , W + !)• (26) 

(1 + e 2 ) 

The results of Np(e) for several values of j3 are shown 
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Figure 1. Differential distribution of stellar orbital eccentricity Np(e) in two extreme cases of mass distribution, such as the point mass 
model on the left panel and the homogeneous model on the right panel. The results are shown by lines for several values of velocity 
anisotropy parameter j3. If Np(e) near e = 1 sensitively changes at some particular value of /3, the results for /3 ± 0.05 are additionally 
shown by dotted lines for the purpose of illustrating its sensitivity. Note that Np(e) is normalized such that L Np(e)de = 1. 



on the right panel of Figure [T] For P < 0.5, Np(e) shows a 
hump-like e-distribution with a single peak at 



-peak 



4(1 - @) - ^13 - 32/3 + 16/3 2 



(27) 



For 0.5 < p < 1 — \/3/4, however, Np(e) has two local 
maxima such as a broad peak at e = epeak and a sharp peak 
at e = 1. Overall behavior monotonically increases with e in 
a range of < e < e pca k, and is kept more or less fiat in the 
range of e pcak < e < 1. For 1 - y/3/4 < P < 1, Np{e) is a 
rapidly increasing function of e. 

For a given value of /3, Np (e) is more weighted at smaller 
e in the homogeneous model, when compared with the point 
mass model. In particular, for P — (isotropic velocity 
dispersion), Np(e) shows a broad hump- like e-distribution 
around a peak at e pca k = 0.36 in the homogeneous model, 
while showing an exactly linear e-distribution in the point 
mass model. This sensitivity, though between two extreme 
cases, could be used to discriminate the likely mass distri- 
bution in more realistic cases to be considered in section 
3. 



2.4 Effect of central mass concentration 

In the cases of central point mass and truncated homoge- 
neous sphere, the shape of np(E,e) is the same as Np(e), 
because np(E, e) is separable in E and e and thus the shape 
of Np (e) is unaffected by g(E) in equation ([4]) . This property 
generally holds when the density distribution in the trun- 
cated sphere is given by p(r) oc 1/r 7 (see Appendix B). The 
homogeneous model in section [2.3.21 corresponds to 7 = 0. 

Using the cases of 7 = 1 (linear potential model) and 
7 = 2 (singular isothermal model) that are intermediate be- 
tween two extreme cases considered above, we can examine 
how Np(e) depends on the central mass concentration. As 
shown in Figure [5] for p — 0, there is a clear trend that the 
e-distribution is peaked at larger e as the halo mass is more 
centrally concentrated. This trend is also true regardless of 




Figure 2. Differential distribution of stellar orbital eccentricity 
Np(e) for four types of model potentials, when the stellar velocity 
dispersion is isotropic {fi = 0). Considered are the homogeneous 
model (7 = 0, section [2.3.21 1. the linear potential model (7 = 1, 
Appendix B.l), the singular isothermal model (7 = 2, Appendix 
B.2), and the point mass model (section l2.3.1ll . in order of increas- 
ing the central mass concentration. There is a clear trend that the 
e-distribution is peaked at larger e as the halo mass is more cen- 
trally concentrated. Note that the distribution is normalized such 
that fg Np(e)de = 1. 



the value of P and is helpful in interpreting the results of 
more realistic models in the next section. 



3 ECCENTRICITY DISTRIBUTION OF HALO 
STARS 

Our formulation in the previous section can apply to more 
general cases of mass distribution, including the isochrone 
model and the NFW model that could form two ends cov- 
ering in-between any realistic cases of mass distribution of 
dark halo. 
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3.1 Energy-dependent eccentricity distribution 
np(e,e) for the isochrone model 

The gravitational potential of the isochrone model (Henon 
1959) is given by 



V(r) 



GM 



b+Vb 2 + r 2 



(28) 



where M is the total mass and b is the scale length parame- 
ter. Obviously, the asymptotic form in the limit of r ^> 6 or 
r<6 approaches the point mass model or the homogeneous 
model, respectively. Thus, this model, though not explaining 
the flat rotation curve of the galaxy disk at greater distances 
from the galaxy center, is important to study the interme- 
diate case of mass distribution by adjusting the scale size of 
the central core. Furthermore, the isochrone model is partic- 
ularly valuable, because fully analytic expression oinp(E, e) 
can be obtained. 

Provided 6 7^ 0, we define useful dimensionless variables 
and effective potential as follows: 



X = 


r 
b' 


e s 


2bE 

~~ GM 


and 








$eff 


(A; 


x) = 


26 
GM 



x = 



2L 2 
bGM ' 



(29) 



nr) + £ 



+ 



A 



1 + Vl + x 2 2x 



Equation ((3| then reads 



ex 2 + 2y/T+x 2 - ( 2 + - ) = 0. 



A 



(30) 



(31) 



This equation has two real and positive solutions if and only 
if 



-1 < e < and < A < A c i r = 



2(1 + ef 



(32) 



We denote the two solutions :r pcr i and a; ap o, and use of them 
gives the relevant quantities in terms of e and e: 



b* 



T - 2 ^GM ( - £r 



L , = bGM _ A 2 



(33) 



(1 + e 4 ) - (1 + e 2 W(l - e 2 ) 2 + 4e 2 e 2 



(34) 



and 



d -£-) E = - b GM(-er 



(1-e 2 ) 



2„2 1 „4 



1 + 2e 2 e 2 + e 



(1 + e 2 ) 



(35) 



>(l-e 2 y + 4e 2 e 2 
Consequently, after tedious algebra, we succeed for the first 



time to obtain analytic expression of flp(e, e) as follows: 



np(e, e) = 87rV6 5 GM(-e)~2 



(1- 



bGM 



1 + 2e 2 e 2 + e 4 



(1-e- 



2\ 2 



(1 + e 2 ) 



- 4e 2 e 2 



-A- 2 - 

e 



+- 



(1 + e 4 ) - (1 + e 2 )J(l - e 2 ) 2 + 4e 2 e 2 



(36) 



We see that n^(e,e) is not separable in e and e. Therefore, 
unlike the point mass and homogeneous models, the shape 
of np (e, e) depends on e as well as /?. Accordingly, derivation 
of Np{e) needs full numerical integration of np(e,e) over e 
with the weight function g(e) specified. 

When f3 = 0, by taking a limit of e, we obtain 



lim nfl=o(e, e) oc e(— e) 2 , 

E— >— 

and 

e hm i n^ ( £ ,e)oc (i + e2)3 



(1+e) 



(37) 



(38) 



These shapes of e-distribution exactly coincide with those 
in equations (I37p and (|38[) . respectively. As understood from 
the definition of e [= 2bE/(GM)], the limit of e — > corre- 
sponds to 6 — ► with E and M fixed, which is equivalent to 
taking a limit to the point mass model. Likewise, the limit 
of e — > —1 corresponds to b — > 00, otherwise such limit of e 
is not attained with E and M fixed, which is equivalent to 
taking a limit to the homogeneous model. 

The shapes of np{e, e) for several values of e and /3 are 
shown in Figure[3] For any value of /3, there is a general trend 
such that eccentric orbits become more and more dominant 
as e increases. However, a marked /3-dependence shows up 
in the shape of np (e, e). 

When /3 < 0.6, np(e,e) has a hump-like e-distribution 
with a peak at e = e pea k. On the other hand, when 0.6 < 
/3 ^ 1, 7i/3 (e, e) has a monotonically increasing e-distribution 
with a peak at e = 1. In particular, for /3 ~ 0.6 and e ~ — 1, 
np(e, e) shows something like a trapezoidal shape, similar to 
the case of ft ~ 0.6 for the homogeneous model (left panel of 
Figure [T}. Furthermore, for /3 > 0.8, highly eccentric orbits 
prominently dominate in the e-distribution. 

In order to understand the situation differently, the 
plots of e pca k at which the e-distribution is peaked for sev- 
eral values of e and j3 are shown on the left panel of Figure [5] 
Here, by taking a limit of e, we can easily confirm, through 
comparison of this figure with Figure [T] that 



lim e 



P cak 



08. e) 



"peak 



03), 



and 



Jin^ e pea k 03, e) = e p °™ k (/3), 



(39) 



(40) 



where superscripts 'pm' and 'horn' correspond to the point 
mass model and the homogeneous model, respectively. More 
generally, when /3 > 0.6, we see that e poa k = 1 for any value 
of e. On the other hand, when /3 ^ 0.5, we see that e poa k is 
an increasing function of both e and j3. 
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Isochrone model p = -1 
2.5 e = -o.9 
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Figure 3. Energy-dependent differential distribution of stellar orbital eccentricity np(e,e) for the isochrone model. In different panels 
for different values of velocity anisotropy parameter 0, shown by lines are the results for dimensionless energy e = —0.9, —0.7, • • ■ , —0.1, 
in steps of 0.2. Note that np(e, e) is normalized such that J„ np{e, e)de = 1. By this normalization, the inclination of n^{e, e) at e = 0, 
which is lower for smaller |e|, helps identify each line. 



3.2 Energy-dependent eccentricity distribution 
n^(e, e) for the NFW model 



Cosmological simulations have been run to reconstruct 
galaxies from the primordial density fluctuations in the uni- 
verse. These numerical results have shown that the dark halo 
has a universal shape of so-called NFW density profile that 
has little dependence on the cosmology (Navarro, Frenk & 



White 1997), such as 
p(r) 



r ) = Po • 



r(a + r) z 



(41) 



where a is the scale length parameter. This density profile 
behaves as p a 1/r for r <§C a, while p oc 1/r for r 3> a. 
The associated gravitational potential is of the form 



V(r) = -4-ivGpoa 



3 ln(l + r/a) 



(42) 
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Provided a / 0, we define dimensionless variables and 
effective potential as follows: 

T 2 

(43) 



E 



and 
$ cff (A;x) 



AiiGpQa 2 



AirGpoa 2 



A 



AiiGpm 4 



V(r) + £ 



ln(l+x) A 
X 2x 2 



(44) 



(45) 



Equation ((3| then reads 

2 A 

ex + x ln(l + x) — — = 0. 

This equation indicates a one-to-one correspondence be- 
tween (e, A) and (e,e), and allows two real and positive so- 
lutions if and only if 



-1 < e < and < A < A cir = x c ln(l + x c ) 



where x c is the solution for 



l + x c 



. 2£= ln(l±£) + , 1 



(46) 



(47) 



1 + a: 

We denote the two solutions a; per i and £ apo (x apo > £ pe ri), 
and use of them gives 



T r = 



1 



2-KGpo 



xdx 



ex 2 + x ln(l + x) 



8ttGpocl [ex t + Xiln(l + Xi)] (xt 



(48) 



£pcri Or Xg^poj^ 

(49) 



and 



dL 2 
de 



-AirGpoa (a; ap o + x pBr i) 



Xpcrf [2e + yq^— J + ln(l + x pcri ) 



p° ( 2e + T+^) + ln ( 1 + Xa P°) 



(50) 



We see that ng(e,e) does not allow analytic expression in 
terms of e and /3. Accordingly, derivation of n^(e, e), as well 
as Np(e) with the weight function g(e), needs full numerical 
integration for the NFW model. 

The results of np(e,e) for several values of e and fi 
are shown in Figure [4] The plots of e pca k at which the 
e-distribution is peaked for several values of e and j3 are 
shown on the right panel of Figure [5] Here, similarly to the 
isochrone model, by taking a limit of e, we can easily confirm 
that 



lim e p8 8k(/3,e) = e^(/3) 
and 



(51) 



t hm i e pcak (/?, e ) = ep p cak (/3), (52) 

where superscripts 'pm' and 'lp' correspond to the point 



mass model and the linear potential model described in Ap- 
pendix B.l, respectively. 

Except for slight shift of the e-distribution to have more 
weight at higher e, overall behavior of n^(e, e) for the NFW 
model is very similar to the isochrone model. Such slight 
shift occurs, because the mass is little more centrally con- 
centrated in the NFW model compared with the isochrone 
model. 

The insensitivity to the choice of gravitational potential, 
as far as it remains realistic, is encouraging, especially when 
our theoretical e-distribution is to be compared with that 
observed for stars in the Milky Way halo. 



3.3 Results of Np(e) 

In the previous subsections, we have derived the in- 
dependent form of np{E,e) for the respective models of 
isochrone and NFW. In order to obtain their eccentricity 
distribution Np(e) in equation ©, we have to specify the 
weight function of g(E) which can in principle be derived 
in a self-consistent way (Lynden-Bell 1962, 1963). Here, in- 
stead of entering into robustness, however, we take a simple 
approximation of g(E) as having the form: 



g(E) = Ae X p[ - — 



(53) 



where A is a constant and a stands for the radial velocity 
dispersion a r ~ 150 km s for the Milky Way halo stars 
(e.g. Yoshii & Saio 1979; Chiba & Beers 2000, 2001). 

We can imagine that halo stars traveling far distantly 
from the galaxy center with near-zero energy would be cap- 
tured by adjacent dark halo. Thus, it is reasonable to in- 
troduce a truncation energy E t above which g(E) should 
vanish. 

The NFW model provides a direct reason to include E t 
in the analysis. The mass of dark halo within the radius r is 
naively given by 



M(r) = 47rp ffl 



hi 



r\ r/a 

a I 1 + r/a 



(54) 



and diverges in the limit of large r. In fact, numerical sim- 
ulations indicate that the NFW density profile applies only 
inside a certain boundary radius but does not apply beyond 
it because of the existence of adjacent dark halos. Such a 
boundary usually used is the virial radius r2oo within which 
the averaged density is equal to 200 times the critical den- 
sity of the universe and the effects by adjacent dark halos 
are negligible. Thus, it is reasonable to place Et at V(r2oo) 
and assume that while halo stars with E < E t stay in the 
system, those with E > E t could be unbound and leave the 
system. 

From all these considerations, we examine how Np{e) 
would be modified with Et taken into account in the anal- 
ysis. Here, we set E t equal to the potential energy V(r2oo) 
and write it in the dimensionless form: 



Et 



ln(l + c) 



(55) 



where c is the concentration parameter defined as c = 
r20o/o.- Use of the kinematic data of the blue horizontal 
branch stars in the Milky Way halo and some CDM simula- 
tions of a halo of M(r2oo) ~ 1O 12 M0 as massive as the Milky 
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Figure 4. Energy-dependent differential distribution of stellar orbital eccentricity ng(e, e) for the NFW model. Others are the same as 
in Figure [3] 



Way halo gives c = 3.9 - 12.5 (Xue et al. 2008), which cor- 
responds to £ t = —0.4 to —0.2. Thus, a choice of this range 
of e t , together with /3 = 0.5 — 0.7 (cf. section 2.2), would be 
appropriate for our analysis of the Milky Way halo. 

We have repeated the calculations of Np (e) for several 
values of a and et in the integration of np(e,e) over e in 
equation ((9J, and find that Np(e) is insensitive to a but sen- 
sitive to Ef The shape of Np(e) is almost the same as that of 
np(et,e). This is because np(et,e) significantly contributes 
to the integration of np(e,e). For example, in a particular 
case of j3 = for the isochrone model, we clearly see such a 



situation from the explicit expression: 



np =0 {e,e)de = 8ttVGM6 5 (1 + ef{-eY 



(56) 



Using the typical combinations of (/3,£t) = (0.5, —0.2), 
(0.5,-0.4), (0.7,-0.2), and (0.7,-0.4) that more or less 
agree with observations of the Milky Way halo, the results 
of Np (e) for both the isochrone and NFW models are shown 
in Figure [S] We see from this figure that as far as reason- 
able values of /3 and et are adopted, the resulting shape of 
Np(e) should be almost linearly proportional to e, except 
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Figure 5. The peak eccentricity e pcak as a function of dimensionless energy e for the isochrone model (left panel) and the NFW model 
(right panel). The results are shown by lines for different values of velocity anisotropy parameter /3. 



for the deviation only at e > 0.7. This is largely regardless 
of adopting either the isochrone model or the NFW model. 
Thus, if the dominant component of the Milky Way halo is 
in dynamical equilibrium, the total eccentricity distribution 
of stellar halo is expected to have a linear trend at e < 0.7 
similar to our results. On the other hand, the behavior of 
predicted iV(e) at e > 0.7, which still shows little differ- 
ence between the isochrone and NFW models, is sensitive 
to /3 and e t - Consequently, such sensitivity can be used for 
a consistency check of the assumed form of f(E,L). These 
predictions in the separate regions of e < 0.7 and e > 0.7 
are testable, given that large kinematical data of halo stars 
are available at present from the SEGUE project or in the 
near future from the Gaia project. 



4 SUMMARY AND DISCUSSION 

Hierarchical clustering scenarios of galaxy formation sug- 
gest that the major merger of at least several subhalos with 
comparable masses would occur at the last stage of galaxy 
formation. This last major merger would cause the violent 
relaxation of halo stars and make them in dynamical equi- 
librium with a dark halo. Based on the assumptions that 
approximate such a status just after the last violent relax- 
ation (section [TJ, we have presented theoretical predictions 
of N(e) for halo stars. This predicted iV(e) should be ob- 
served for the Milky Way halo if it is an isolated system and 
the subsequent variation of the potential is quiescent enough 
to conserve the eccentricity of each star. 

However, recent nearby observations suggest that at 
least some part of the Milky Way halo may have origi- 
nated from accreted satellites, which possibly deviates the 
observed iV(e) from our predictions. For example, if infalling 
satellites break up and spread their member stars into the 
field, these stars would show peculiar eccentricity distribu- 
tion which necessarily imprints the initial condition of the 
progenitor satellites. In addition, if such satellites locally 
disturb the halo potential, some in-situ halo stars may have 
altered their orbits (e.g. Zolotov et al. 2009). With an in- 
vention of segregating in-situ halo stars from infalling stars, 



we might be able to well understand the nature of accretion 
and distortion of satellites. 

Numerous authors subdivided halo stars into some 
'components' and examined the correlations between chem- 
istry, age and kinematics of stars in each component. Carollo 
et al. (2010) obtained reliable eccentricities for ~ 10, 000 
halo stars within 4 kpc of the sun and decomposed them 
into the inner and outer halo components having distinct 
eccentricity distributions from each other. Since their sam- 
ple is local and is inherently biased in favour of stars that 
stay longer in the surveyed region, our formalism, which is 
designed to predict N(e) of the whole stellar halo, has to be 
modified for the purpose of fair comparison with their data. 
Through proper incorporation of effects of such a bias, we 
can still predict N(e) for a local sample by fully taking into 
account a probability of finding each of halo stars in the 
surveyed region. This will be done in a separate paper in 
preparation. On the other hand, our formalism can directly 
apply to a global, and therefore less biased, sample of halo 
stars with reliable orbital eccentricities, such as those from 
next generation surveys including the Gaia mission. In either 
case, the analytical approach in the present paper certainly 
forms a basis that serves as a useful tool for analysing the 
kinematics of the stellar halo. 

Large, unbiased database of halo stars would enable us 
to test whether a given component is in dynamical equi- 
librium by comparing the observed and predicted shape of 
N(e). Such comparison would hopefully discover some re- 
laxed components, and their adiabatically conserved shape 
of N(e) would carry some useful information of the physics 
of violent relaxation. Moreover, the spatial distribution of 
these relaxed components would enable us to see how far 
out in the halo the violent relaxation has exerted and how 
strongly it has affected the stellar halo. If the information 
of last violent relaxation, yet to be known observationally, 
is gained in this way, more precise assessment to the early 
evolution of the Milky Way would be possible, and our un- 
derstanding of its formation would greatly be advanced. 

Our current calculations of N(e) are certainly very sim- 
ple and can be improved by using more realistic assump- 
tions. For example, we can modify our analysis to allow 
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Figure 6. Differential distribution of stellar orbital eccentricity Np(e) for the isochrone model (left panel) and the NFW model (right 
panel). Adopted combinations of (/3, et) more or less agree with observations of the Milky Way halo. Note that the normalization factor 
of Np(e) is arbitrarily chosen so that the nearly linear trend op to e rj 0.7 is clearly seen. 




axisymmetric potentials including a disk-like component as 
well as a bulge. Preliminary analysis has confirmed that 
inclusion of a disk-like component would cause no signifi- 
cant change in the linear trend of iV(e) described in section 
I3.3I which will be discussed in a separate paper. Also, our 
choice of f(E, L) having the form in equation Q has to be 
extended to allow the radial dependence of j3(r). Further 
elaborate modeling of N(e) with these theoretical improve- 
ments, when applied to future large survey of halo stars, 
would then provide a promising way of unraveling myster- 
ies of the galaxy formation and evolution in a paradigm of 
hierarchical clustering in the ACDM cosmology. 
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APPENDIX A: ALLOWED REGION OF (E,L) 
FOR BOUND ORBIT 

A steady, bound orbit in a gravitational potential V(r) 
generated by a density distribution p(r) is only possible 
in a subset of energy E and angular momentum L that 
allows two real and positive solutions for equation ((3|. We 
discuss such an allowed region of (E, L) in this appendix. 
We begin with the effective potential 



V aa (L;r) = V(r) + ^. 

Then, from the definition, we obtain 



J^ ff (L;r))^l[GAf(r)r-L 2 ] 



(Al) 



(A2) 
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where M(r) is the total mass inside the radius r. Since 
GM(r)r is a monotonically increasing function of r and 
it satisfies 

lim [GM(r)r] = 0, and lim [GM(r)r] = oo, (A3) 

r— >0 7' — ^oc 

there always exists an radius r c = r c (L) for which 
GM(r c )r c — L 2 , thus yielding 



^V aS {L;r)\ ^0, ifr^r c 



Since we have 
d 



d(L 2 
and 

d 



r c = [G (47rr^(r c ) + M{r c ))] 1 > 0, 



1 



(A4) 



(A5) 



(A6) 



d(L2) K ff (L;. c (L)) = ^>0, 

the allowed range of E with L fixed can be expressed as 

V cB (L;r c {L))<E<0. (A7) 

Here, we define the zero of V(r) so that lim, -^oo V(r) = 0. 
Thus, for any given L, we obtain 



lim V eB (L;r) = 0, 



(A8) 



which validates that the upper bound of inequality ()A7|) 
should be zero. As for the allowed region of L when E is 
fixed, we obtain 



< L< L cil {E), 

for which Ldr(E) is the solution of 

E = V e s(Lcir;r c (L clT )). 



(A9) 
(A10) 



and we will refer to this potential a 'truncated linear po- 
tential.' We consider only stars with E < E t = —GM/rt, 
which guarantees the stars to be confined inside the trun- 
cated radius r t . Thus, bound orbits within the truncated 
sphere are allowed if -E m in < E < E t where we note 
E m in = 2E t . In this limited range of E, there are two real 
and positive solutions for equation (3), or equivalently, 

2GMrJy\ - 2r 2 (E - E min ) (I J + L 2 = 0, (B4) 



if and only if 
< D < 2, 

where 

27G 2 M 2 L 2 



<±ri{E - E min ) 



(B5) 



(B6) 



In this allowed region, two real and positive solutions for 
equation (|B4|) are as follows: 



2r 
n = „\, (E - E min ) Xi (i = apo or peri; r apo > r pcri ), 



with Xapo and :r pe ri given, respectively, by 



a; ap o = — + cos#, and a; pcr i = — + cos 



thus 



47T 



(B7) 



(B8) 



APPENDIX B: OTHER MODELS OF 
TRUNCATED MASS DISTRIBUTION 

We present the derivation of Np (e) for two models of trun- 
cated power-law mass distribution: 



(3- 7 )M I r 
p(r) = < ^ r t V r » 

o, 



(7 < 3) if r < r t 
otherwise, 



(Bl) 



cos 6- cos [^- + 0\ 
1 + cos 9 + cos [fe + 1 



where 



Itan-(^ 



tan 



\/2D-D 2 



if < D < 1 



if f < D < 2. 



(B9) 



(B10) 



where M is the total mass of the dark halo and r t is the 
truncation radius. We note that the truncated homoge- 
neous model presented in section 12.3.21 is a special case of 
7 = in equation (|B1[) . 

Bl Linear potential model (7 = I) 

The NFW density profile has a central cusp and behaves 
like p(r) oc 1/r in the limit of Small r. This density profile 
corresponds to 7 = 1 in equation (|B1|I . and we have 



p(r) 



M 

2tt 



V?\n) 



0. 



if r < n 

otherwise. 



(B2) 



The gravitational potential arising from this density pro- 
file is given by 



V(r) = 



2GM 

rt 
GM 



+ ^j, ifr<r t 

otherwise, 



(B3) 



Consequently, D has a one-to-one correspondence to e, so 
with 0, Xapo, and a; por i. Use of these quantities gives 



27G 2 M 2( minj ' 



T r = 2\/3 



rj ^E - E min 
GM 



xdx 



i- 



r3 4- 3^-2 _ D 



(BH) 



= , (BI2) 



and 
dL 2 \ 4r t V2D - D 2 



de J E 9G 2 M 2 



(E — -Emin) 



(f + cos6» + cos [±f 



sin 9 (I + 2 cos [*f + 6] ) - sin [4f + 0] (I + 2 cos ( 



(B13) 
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By these expressions, we obtain 

^ „ r 2 / Ar 4 



GM\21G 2 M 2 1 
D~ l3 y / 2D~D 2 (E - E u ^)i~ z(i 

sin 6 (1 + 2 cos [^ + 9] ) - sin [^ + 0] (1 + 2 cos 0) 



xdx 



_t3 i 3 2 _ D 



(B14) 



Since 6, D, £ apo , and x pcr i depend only on e, np(E,e) is 
separable in _E and e, so that 



iV>(e) = 24v / 3tt 2 ^- f— i^- 



4 x J ' 

GM \ 27G 2 M 2 



D~ I3 < S /2D~D 2 



g(E)(E-E min )*-^dE 



(1 + cos 6 + cos [^ +1 



sin # ( 1 + 2 cos 



+ i 



sin [^-+0] (l + 2cos0) 
:rd:r 



D_ 

4 



(B15) 



Therefore, the shape of Np{e) is not affected by E or 
g(E), like the point mass model and the truncated model 
with any 7. The results of Np(e) in the linear potential 
model are shown on the left panel of Figure IB1I We see 
that Np (e) is a monotonically increasing e-distribution for 
0.52 < /3 < 1, and Np(e) has a hump-like e-distribution 
with a single peak for /3 < 0.5. In particular, Np = o(e) 
reaches its maximum at e pca k = 0.43. In the intermediate 
range of 0.5 < /3 < 0.52, Np(e) shows something like a 
trapezoidal e-distribution, which shows a monotonically 
increasing e-distribution for < e < e pca k and a more or 
less flat behavior for e pca k < e < 1, where e pca k — 0.8. 



are allowed if — 00 < E < E t . In this range of E, there 
are two real and positive solutions for equation ([3]), or 
equivalently, 



-2iLYr-K,'.Ur, )(£-j + 2GMrJ^\ In (y ] -:- L~ - 0. 



if and only if 



0<D < 



2exp(l)' 



where 
D 



2GMn exp [2 (l + gf ) 



(B18) 



(B19) 



(B20) 



In this allowed region, two real and positive solutions are 
as follows: 



r - n «'xi ) ( 1 + -£— ) Bi, (i = apo or peri; r apo > r pc 



with a; apo and x por i are the solutions for 



x lnx + D = 0. 



(B21) 



(B22) 



By this equation, D has a one-to-one correspondence to 
e, so with x P eri and a; apo . Use of these quantities gives 



T r 



2± 
GM 



exp 



GM 1 



xdx 



, ri y/—D — x 2 lnx 



(B23) 



B2 Singular isothermal model (7 = 2) 

One of the most strong constraints on the gravitational 
potential of the halo is that it has to be consistent with 
the observed flat rotation curve of galaxy disk. In this 
sense, the truncated singular isothermal model, which au- 
tomatically reproduces the flat rotation curve in the radial 
range of < r < r t , is said to be one of the simple and 
realistic models. The density profile of this model is given 
by 



p{r) = 



0, otherwise, 



M 

4tt 



(B16) 



which corresponds to 7 = 2 in equation (|B1|) . The gravi- 
tational potential arising from this density profile is given 
by 



V(r) = 



GM 



if r < rt 

otherwise. 



(B17) 



We consider only stars with E < E t = —GM/r t , which 
guarantees the stars to be confined inside the truncated 
radius r t . Thus, bound orbits within the truncated sphere 



L" = 2GMr t Dexp 

and 
dL 2 



1+^1 

GM 



de 



= —GMrt exp 



1 + ^ 



GM J 



l^apo -\- XpcriJ 



2D x 2 



2D 



Consequently, we obtain 



n (i (E,e) = 4V2-K 2 \j GMrl[2GMr t D]~ 
x exp 



^-^(l + g§) 



X 



(Xapo ~T~ XpcriJ 



2D 4 cli - 2D 
Xapo xdx 



ori V-D-x 2 lm 



(B24) 



(B25) 



(B26) 
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Figure Bl. Differential distribution of stellar orbital eccentricity ATg(e) in two cases of truncated mass distribution, such as the linear 
potential model (7 = 1) on the left panel and the singular isothermal model (7 = 2) on the right panel. The results are shown by lines 
for several values of velocity anisotropy parameter /3. If Np(e) near e = 1 sensitively changes at some particular value of 0, the results 
for /3 ± 0.05 are additionally shown by dotted lines for the purpose of illustrating its sensitivity. Note that Np (e) is normalized such that 

/ 1 iV j9 (e) + «fe = l. 



Since D, x pcT i, and x^ po depend only on e, np{E,e) is 



separable in E and e, so that 



Np[e) = 4\/2tt 2 ^GMrl [2GMr t D] ' 



g(E)exp 



(3 - 2/3) 1 + 



nE \ 
GMJ 



dE 



\3^apo \ ^pcrij 



2D 



-2D 

xdx 



y/—D — x 2 \nx 



(B27) 



Thus, the shape of Np(e) is not affected by g(E), like 
the point mass model and the truncated model with any 
7. The results of Np{e) in the singular isothermal model 
are shown on the right panel of Figure IB1I We see that 
Np (e) shows a monotonically increasing e-distribution for 
/3 > 0.45, while having a single peak for /? < 0.45. 



